UNCLASSIFIED 


_ AD  NUMBER _ 

AD826609 

LIMITATION  CHANGES 
TO: 

Approved  for  public  release;  distribution  is 
unlimited. 


FROM: 

Distribution  authorized  to  U.S.  Gov't,  agencies 
and  their  contractors ; 

Administrative/Operational  Use;  11  JAN  1968. 
Other  requests  shall  be  referred  to  Air  Force 
Technical  Applications  Center ,  Washington ,  DC 
20333. 


_ AUTHORITY 

AFTAC  per  DTIC  form  55 


THIS  PAGE  IS  UNCLASSIFIED 


rRUMENTS  INCORPOR 
ice  Services  Division 
P.  O.  Box  5621 
las,  Texas  75222 


Contract  No.  AF  33*657) -16678 


LARGE -ARRAY  SIGNAL  AND  NOISE  ANALYSIS 


Quarterly  Report  No.  6 
1  October  1967  through  31  December  1967 


Frank  H.  Binder,  Program  Manager 


TEXAS  INSTRUMENTS  INCORPORATED 
Science  Services  Division 
P.  O.  Box  5621 
Dallas,  Texas  75222 


Contract  No.  AF  33(657)-l6678 


Prepared  for 


AIR  FORCE  TECHNICAL  APPLICATIONS  CENTER 
Washington,  D.  C.  20333 

Sponsored  by 


ADVANCED  RESEARCH  PROJECTS  AGENCY 
ARPA  Order  No.  599 
AFTAC  Project  No.  VT/6707 


1 1  January  1 968 


acioncs  aarvicas  division 


Texas  Instruments 

INCORPORATED 
SCIENCE  SERVICES  DIVISION 

1 1  January  1968 


Air  Force  Technical  Applications  Center 
VELA  Seismological  Center 
Headquarters,  USAF 
Washington,  L>  C.  20333 


Attention:  Captain  Carroll  F.  Lam 

Subject:  Sixth  Quarterly  Report  Covering  Period  October  1,  1967 

Through  December  3 1 ,  1967. 


AFTAC  Project  No.  : 
Project  Title: 

ARPA  Order  No.  : 

Name  of  Contractor: 

Date  of  Contract: 

Amount  of  Contract: 
Contract  Number: 
Contract  Expiration  Date: 
Program  Manager: 


VT/6707 

Large  Array  Signal  and  Noise  Analysis 

599 

Texas  Instruments  Incorporated 

16  May  1966 

$1, 083, 696 

AF  33(657)-  16678 

25  June  1968 

Frank  H.  Binder 

Area  Code  214 

238-3473 


Gentlemen: 

Approval  of  the  one  year  extension  in  the  amount  of  $490,  023 
was  concluded  during  the  last  quarter. 

Below  is  set  forth  the  work  progress  against  the  major  tasks 
remaining  under  the  contract. 

PUBLICATION  OF  SPECIAL  REPORTS  DEALING  WITH  WORK  PREVIOUSLY 
COMPLETED 

The  following  special  reports  dealing  with  previously  completed 
work  were  published  during  the  past  quarter. 
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LASA  Special  Report  No.  3  -  SUBARRAY  PROCESSING. 

LASA  Special  Report  No,  6  -  ANALYSIS  OF  SUBARRAY 
WAVENUMBER  SPECTRA. 

LASA  Special  Report  No.  12- ANALYSIS  OF  LONG  PERIOD 
NOISE. 

LASA  Special  Report  No.  13-NOISE  COHERENCE  AMONG 
SUBARRAYS. 

LASA  Special  Report  No.  14-WIENER  NON-TIME 
STATIONARY  PROCESSING. 

The  following  reports  are  in  different  stages  of  publication  and 
should  be  completed  by  February,  1968. 

LASA  Special  Report  No.  4  -  SPACE  AND  TIME  VARIABILITY 
OF  THE  LASA  NOISE  FIELD. 

LASA  Special  Report  No.  10  -EQUALIZATION  STUDIES. 

LASA  Special  Report  No.  1 1 -RESOLUTION  OF  EVENTS. 

SHORT  PERIOD  NOISE  ANALYSIS 

K-line  spectra  have  been  run  on  two  summer  noise  samples.  The 
spectra  have  not  yet  been  analyzed.  These  spectra  should  give  an  indication  of 
any  fundamental  seasonal  changes. 

A  report  on  the  summer  noise  analysis  will  be  published  during 
the  coming  quarter.  This  report  will  describe  the  wavenumber  spectra  of  the 
summer  noise  samples  and  compare  them  with  the  previously  published  results 
for  the  winter  noise  samples.  (LASA  Special  Scientific  Report  No.  6,  ANALYSIS 
OF  SUBARRAY  WAVENUMBER  SPECTRA.) 

A  study  of  locally  generated  noise  at  the  LASA  will  be  begun 
in  February.  The  purpose  of  this  sUdy  will  be  to  try  to  identify  and  describe 

noise  sources  within  the  large  array. 
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LONG  PERIOD  NOISE  STUDIES 

•  Analysis  of  Summer  Noise  Samples 

The  purpose  of  this  study  is  to  compare  the 
characteristics  of  the  long-period  noise  field 
during  the  summer  months  with  those  observed 
during  winter  months  (Special  Report  No.  12). 

Five  noise  samples  will  be  analyzed;  these  have 
been  demulitplexed,  transferred  to  the  System 
360,  and  are  presently  being  despiked  and 
plotted.  Analysis  techniques  will  parallel  those 
described  in  Special  Report  No.  12. 

•  Further  Study  of  the  Noise  Below  0.  05  cps 

The  long  noise  sample  to  be  used  to  study  the  noise 
below  0.  05  cps  consists  of  five  back-to-back 
80 -minute  tapes.  These  tapes  have  been  de¬ 
multiplexed,  transferred  to  the  System  360,  and 
are  presently  being  despiked  and  plotted.  Plots 
have  shown  that  the  first  tape  has  a  small  sur¬ 
face  wave  arrival  (primarily  on  the  verticals), 
and  the  fifth  tape  has  an  89-point  segment  of  bad 
data.  The  feasibility  of  including  part  or  all  of 
these  tapes  in  the  long  noise  sample  is  presently 
being  investigated.  This  long  noise  sample  will 
allow  more  careful  estimates  of  the  low  frequency 
spectra  and  coherence. 

In  addition,  a  program  is  being  written  to  demulti¬ 
plex  the  microbarographic  data  associated  with  the 
long  noise  sample.  It  is  anticipated  that  preliminary 
data  handling  (including  the  microbarographic  in¬ 
formation)  will  be  completed  early  in  February, 

1968. 

The  microbarograph  data  will  be  used  to  check 
for  correlation  between  pressure  fluctuations  and 
long  period  seismometer  output  at  low  frequencies 
(f  <  .  05  cps). 
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Continued  Studies  of  Love-Wave  Noise  Ene rgy 

Past  studies  (LASA  Special  Scientific  Report 
No.  12)  have  indicated  that  point-like  sources 
of  Rayleigh  mode  energy  may  also  be  sources 
of  Love  mode  energy. 

Some  horizontal  components  have  been  rotated 
in-line  and  transverse  to  the  very  strong  storm 
source  indicated  in  the  February  7  noise  sample. 
Power  spectra  and  co1  fences  have  been  generat¬ 
ed  to  see  if  there  are  indications  that  this  storm 
source  is  generating  Love  mode  energy. 

In  the  region  near  the  7-sec  microseism  peak, 
the  horizontal  components  are  generally  5-10  db 
noisier  than  the  vertical  components  (LASA  No. 

12).  This  may  be  the  result  of  (1)  Love  mode 
energy,  (2)  non-seismic  noise,  (3)  change  in 
excitation  of  horizor.tal  to  vertical  motion  of 
Rayleigh  mode  energy. 

Possibility  No.  2  is  unlikely  based  on  the 
apparent  coherence  of  the  horizontal  components 
in  this  frequency  range  (LASA  No.  12).  Possi¬ 
bility  No.  3  can  be  fairly  well  confirmed  or 
denied  by  calculating  the  excitation  functions  for 
the  LASA  crust.  These  calculations  are  being 
made  and  will  be  completed  and  reported  during 
the  coming  quarter. 

Study  of  Vertical  Component  Noise  Not  Predictable 
From  Horizontal  Components 

In  general  about  90%  of  the  noise  power  on  the 
verticals  is  predictable  from  nearby  horizontals 
at  the  15-18-sec  microseism  peak.  In  order 
to  try  to  dissect  the  remainder  of  the  vertical 
component  power,  a  simultaneous  prediction 
of  all  verticals  from  all  horizontals  will  be  done. 

This  operation  is  conveniently  accomplished  in 
the  frequency  domain  from  the  crosspower  matrix 
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of  all  channels  (vertical  and  horizontal). 

This  operation  yields  a  conditional  covariance 
only;  that  is,  the  covariance  matrix  of  the 
unpredictable  (from  horizontals)  part  of  the 
vertical  component  noise.  The  resulting 
covariance  matrix  can  be  used  for  coherency 
and  wavenumber  analysis. 

Using  such  a  large  number  of  channels 
(potentially  21V,  21E  and  21N)  requires  a  very 
long  noise  sample  to  obtain  adequate  statistics. 
The  long  noise  sample  described  above  will  be 
used  for  this  study.  Depending  on  how  long  a 
sample  is  finally  usable,  the  number  of  channels 
will  be  adjusted  to  give  adequate  statistics. 

Stability  of  Long  Period  Noise 

Time  domain  signal  extraction  and  prediction 
filters  were  designed  on  one  noise  sample  and 
applied  to  the  design  interval  and  to  another  time 
gate  separated  by  about  6  hours. 

This  data  is  being  combined  with  a  study  of  the 
usefulness  of  horizontals  in  P-wave  extraction 
which  is  not  yet  complete.  A  report  on  this 
study  will  be  forthcoming  during  February. 

Correlation  of  Storms  at  Sea  with  the  Long 
Period  Noise  Sources 


A  study  of  the  correlation  between  the  long 
period  noise  sources  (as  seen  in  wavenumber 
spectra)  and  sea  conditions  has  been  completed. 
The  correlation  appears  to  be  meaningful  but 
many  details  are  perplexing.  The  results  are 
detailed  in  LASA  Special  Scientific  Report  No. 

17  which  has  been  submitted  to  AFTAC  for 
approval. 
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LONG  PERIOD  SIGNAL  ANALYSIS  AND  DISSECTION 

A.  EXTRACTION  OF  SIGNALS  FROM  AMBIENT  NOISE 

The  goal  of  this  work  is  to  evaluate  various  techniques  for 
extracting  the  bodywave  phases  and  surface  wave  nodes  of  an  event  in  the 
presence  of  ambient  noise.  Both  off- and  on-line  procedures  will  be  studied. 

In  the  off-line  case  we  assume  that  the  approximate  epicenter  and  P-wave 
arrival  time  are  known.  For  this  case,  a  number  of  specific  questions  arise. 

1)  What  is  the  beam  steer  performance? 

2)  Can  MCF  provide  superior  performance  or  equivalent  performance 
with  reduced  array  size  and  number  of  elements? 

3)  How  does  microseismic  storm  noise  affect  our  results? 

4)  How  useful  are  the  horizontals? 

To  gain  some  insight  into  the  problem,  70-minute  noise  samples 
from  December  2  and  3,  1966  have  been  analyzed.  First  the  following  averages 
(straight  sums)  of  vertical  traces  were  computed. 

9  channels  -  AO,  C  ring  and  D  ring 
12  channels  -  AO,  C  ring,  D  ring,  and  E  ring  except  El 

(Dec.  2)  18  channels  -  all  except  El,  FI,  and  F3 

(Dec.  3)  20  channels  -  all  except  FI 

Next,  infinite  velocity  signal  extraction  filters  were  designed 
for  the  first  two  combinations  listed  above.  The  December  3  noise  sample  was 
used  to  develop  noise  statistics  for  the  filter  design,  and  the  filters  were  applied 
to  both  noise  samples.  Ratios  of  the  AO  noise  spectrum  to  the  spectra  of  the 
outputs  of  the  various  processors  were  computed.  The  results  will  be  published 
in  a  special  report.  The  following  points  of  interest  are  noted. 

1)  Noise  suppression  by  the  beam  steer  processors,  while  differing 
in  small  details,  is  essentially  the  same  for  these  two  noise 
samples. 

2)  In  the  range  0.05  to  0.20Hz  for  the  December  3  sample,  each 
MCF  is  everywhere  superior  to  the  corresponding  beam  steer. 
Either  MCF  is  superior  to  any  of  the  beam  steer  processors 
except  for  a  small  range  near  0.  14  Hz  where  the  20 -channel 
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average  outperforms  them. 

3)  Performance  of  the  MCF  processors  on  the  December  2 
samplp  '.s  about  2  -  3  db  poorer  than  on  the  December  3 
sample. 

In  addition,  a  prediction  error  MCF  was  designed  using  the  AO,  Bl,  B3, 
and  B4  horizontals  to  predict  the  AO  vertical.  This  processor  was  notably 
inferior  to  any  of  the  others  in  suppressing  the  noise.  The  reasons  for  this 
are  not  fully  understood. 

Finally  an  infinite  velocity  signal  extraction  MCF  was  designed 
using  all  the  AO  and  B  ring  sensors  except  the  Bl  vertical  and  B2  horizontals. 

The  result  was  just  slightly  inferior  to  the  9  and  12  channel  MCF  systems. 

Since  the  B  ring  horizontals  were  not  effective  in  predicting  the  AO  vertical 
this  suggests  using  the  AO  and  B  ring  verticals.  This  is  in  progress  but  has 
not  been  completed. 

These  results  indicate  that  multi-channel  filtering  of  the  inner 
rings  may  be  as  effective  as  beam  steering  the  entire  array  for  bodywave  signal 
extraction.  This  will  be  further  investigated  by  actually  performing  the  extraction 
on  recorded  signals.  In  designing  the  MCF  processors  there  is  some  question 
regarding  the  best  way  to  generate  the  noise  statistics.  Ordinarily,  it  would 
seem  best  to  use  a  noise  sample  immediately  preceding  the  event.  If,  however, 
there  is  overlap  between  various  phases  of  the  event,  the  use  of  a  time  gate 
including  the  entire  signal  might  result  in  better  extraction  of  the  individual 
phases.  This  problem  will  be  studied. 

The  type  of  analysis  performed  on  the  2  December  noise  samples 
will  be  repeated  as  an  initial  approach  to  the  problem  of  surface  wave  extraction. 
Here  rather  than  straight  sum  we  will  use  time  shift  and  sum,  and  the  MCF 
processors  will  be  designed  for  c  surface  wave  signal  model.  In  this  way  we 
hope  to  get  an  indication  of  which  types  of  processing  are  apt  to  be  effective  on  this 
problem.  The  December  noise  samples  used  in  the  above  studies  were  domina¬ 
ted  by  point-like  noise  sources  indicating  that  they  contained  a  substantial  amount 
of  microseismic  storm  noise.  The  presence  of  such  noise  is  probably  more 
significant  in  the  surface  wave  extraction  problem  than  for  bodywave  extraction. 
Therefore,  this  analysis  will  include  noise  samples  representing  the  quiet 
ambient  noise  case  and  the  storm  noise  case. 
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B.  EXTRACTION  OF  SIGNAL  FROM  STRONG  INTERFERING  EVENT 

This  effort  deals  with  the  extraction  of  event  phases  which 
arrive  simultaneously  with  a  strong  signal  from  another  event.  The  problem 
is  to  suppress  the  interfering  signal  while  retaining  and  dissecting  the  target 
event.  Several  off-line  approaches  to  the  problem  will  be  studied. 

1)  Beam  steer.  This  technique  is  straightforward.  In  the  case 
of  the  surface  wave  phases  of  the  target  event,  consideration 
of  dispersion  in  doing  the  beam  steer  should  provide  some 
improvement  over  simple  time  shift  and  sum.  This  effect, 
however,  is  expected  to  be  small. 

2)  Another  method  is  to  time  shift  all  n  traces  to  make  the  inter¬ 
fering  signal  align  (again  possibly  taking  dispersion  into 
account).  Then  by  subtracting  one  trace  from  all  the  rest  one 
obtains  n-1  traces  containing  distorted  target  signal  and  ambient 
noise.  The  signal  is  reconstructed  and  ambient  noise  reduced 
by  subsequent  MCF. 

If  the  signal  and  the  interfering  signal  exactly  fit  a  plane  wave 
hypothesis  one  could  entirely  eliminate  the  interfering  signal 
while  passing  the  desired  signal.  Obtaining  very  large  rejection 
probably  depends  on  forcing  the  best  possible  fit  to  a  plane  wave 
model.  Interpolation  for  time  shifts  of  less  than  one  second  and 
perhaps  frequency  dependent  equalization  will  be  attempted. 

3)  The  third  method  is  to  design  a  MCF  to  operate  directly  on  the 
seismic  traces.  A  problem  here  is  to  determine  the  best  way 
to  generate  noise  statistics  for  the  filter  design.  One  way  is 
to  use  a  theoretical  noise  model.  If  one  uses  measured  noise 
the  question  is  whether  tv  use  a  time  gate  containing  all  of  the 
interfering  event  or  just  that  portion  which  precedes  the  arrival 
of  the  target  event. 

4)  A  fourth  technique  which  may  be  useful  is  to  form  one  beam 
toward  the  interfering  signal.  One  then  uses  the  latter  to  pre¬ 
dict  the  former  and  calls  the  prediction  error  the  signal 
estimate.  A  technique  of  this  type  would  probably  be  easy  to 
imple  lent  on-line. 

In  view  cf  the  nature  of  the  principal  noise  component  in  this 
problem,  it  appears  likely  that  results  superior  to  those  realized  by  a  beam  steer 
processor  should  be  available.  In  order  to  test  this  hypothesis  it  is  necessary 
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to  first  do  the  beam  steer  and  then  to  perform  other  types  of  processing. 

To  accomplish  this  the  Rayleigh  phase  from  a  California  event  has  been 
imbedded  in  the  Rayleigh  phase  from  a  Greenland  Sea  event  with  its  peak 
amplitude  scaled  to  half  that  of  the  Greenland  Sea  event.  Epicentral 
directions  were  225°  and  20°  respectively.  A  simpie  time  shift  and  sum 
aimed  at  the  California  event  resulted  in  about  a  4:1  improvement  in  the 
amplitude  of  that  event  over  the  interfering  one.  This  work  will  be  continued 
by  evaluating  more  complex  approaches  to  the  problem.  Similar  analyses 
will  be  performed  by  compositing  events  with  lesser  amounts  of  azimuthal 
separation. 

C.  OTHER  STUDIES 

The  lateral  refraction  of  surface  waves  is  a  matter  of  con¬ 
sequence  in  this  work.  This  results  in  propagation  vectors  at  the  array  which 
differ  from  the  great  circle  path  from  the  epicenter  to  the  array.  If  effective 
enhancement  of  the  signal  is  to  be  accomplished  the  propagation  vectors  for 
the  signal  and  any  interfering  events  must  be  determined.  A  technique  employ¬ 
ing  the  outputs  of  three  seismometers  to  estimate  the  speed  and  direction  as  a 
function  of  frequency  has  been  developed.  This  has  been  applied  to  the 
California  and  Greenland  Sea  events.  Analysis  of  other  events  in  our  library 
will  be  performed. 

The  study  of  high  resolution  wavenumber  spectra  for  the  separa¬ 
tion  of  event  phases  has  been  continued.  A3  reported  previously  this  technique 
has  been  applied  to  an  event  which  has  clear,  time-separated  compressional, 
shear,  and  Rayleigh  phases.  When  the  time  gate  used  for  computing  the  spectra 
included  the  Rayleigh  phase,  the  spectra  did  not  show  the  shear  energy.  Sub¬ 
sequently  this  has  been  repeated  with  a  data  gate  containing  only  the  shear  phase. 
In  this  case,  the  shear  energy  was  well  defined  in  the  resulting  f-Tc  plots. 
Continued  work  in  this  area  is  planned. 

HIGH  RESOLUTION  WAVENUMBER  SPECTRA  RESEARCH 


Probabilistic  processing  in  the  frequency-wavenumber  domain  is 
developed.  Conventional  high  resolution  frequency-wavenumber  spectra  are 
shown  to  be  equivalent  to  a  probabilistic  processor  when  the  noise  field  is  assumed 
to  be  Gaussian  and  uncorrelated  between  sensors.  Coherent  noise  fields  are 
shown  to  reqx  ire  more  complex  processors  for  best  results. 

A.  DERIVATION  OF  THE  PROBABILISTIC  PROCESSING  EQUATION 

Frequency  domain  probabilistic  processing  is  a  method  \ising  seismic 
array  data  to  detect  plane  wave  signals  in  the  presence  of  ambient  seismic  noise. 
This  method  is  based  on  the  assumption  that  the  array  data  is  Gaussian  with 
mean  zero  and  known  crosspower  matrix  CQn3  or  CQoVnI  . 
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depending  on  the  absence  or  presence  of  signal.  The  array  output  is  con¬ 
sidered  to  be  K  complex  numbers  resulting  from  the  collateral  Fourier 
transform  of  time  gates  from  K  channels  of  simultaneous  array  time  domain 
data. 


The  criteria  for  estimating  the  presence  or  the  absence  of  a 
plane  wave  signal  is  analogous  to  testing  the  simple  hypothesis  that  the 
observed  data  has  a  crosspower  matrix  t QhI  with  the  alternate  multiple 
hypotheses  that  the  data  is  from  a  set  of  crosspower  matrices  of  the  form 
[Q's'+n!  •  The  superscript  m  refers  to  a  particular  plane  wave  signal  model. 

The  derivation  of  the  Gaussian  multichannel  sampled  data  pro¬ 
cessing  equation  staru  with  the  multi-variate  Gaussian  probability  distri¬ 
bution.  1 


p (*.,**,. [XI  Mi}. 

In  this  equation,  matrix  notation  is  used  as  follows: 

ttlt  or  X*.  =4,  set  of  complex  data  values  Xi j  Y-t  j  •  ••  /  %  k 

with  Xx  being  the  Fourier  transform  of  the 
j time  gate. 

IX]  or  X  =  crosspower  matrix  of  the  data,  i.e.  , 


The  superscript  T  designates  the  matrix  transpose;  the  superscript  -/ 
designates  the  inverse  matrix;  and  the  superscript  *  signifies  the  matrix 
complex  conjugate.  Vertical  I  ines  bordering  a  square  matrix  indicates  the 
matrix  determinant. 

In  Equation  II- 1,  K  is  the  number  of  data  values.  In  a  multi¬ 
channel  problem,  K.  would  usually  be  equal  to  the  number  of  channels. 

The  detection  of  a  plane  w.ve  signal  in  the  presence  of  ambient 
seismic  noise  can  be  formulated  as  a  multi-possibility  multichannel  Gaussian 
problem.  That  is,  one  can  assume  that  the  data,  Xi.  is  a  sample  of  a  stationary 
Gaussian  time  series  whose  K  by  K  crosspower  matrix  is  noise,  CQmI  ,  or  one  of 
a  set  of  signal  plus  noise  matrices,  [QC+mI  • 
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To  completely  specify  the  problem  using  a  Bayesian  approach, 
the  a  priori  probabilities  PCQn),  PCQ«?  +  vO>  of  the  possibilities  must 
be  known  or  estimated.  The  problem  is  then  to  obtain  the  a  posteriori 
probabilities  of  noise  alone,  P  (  Qw  i  \Ci),  and  of  signal  plus  noise, 

WQ*4NlXi^by  analyzing  the  data  X*.  It  is  understood  that  CQ«v,^l  is  to  range 
over  all  allowable  plane  wave  signal  models. 


With  regard  to  the  seismic  detection  problem,  Q*j  stands  for 
ambient  seismic ^noise  and  P(Qn1  would  be  the  fraction  of  time  that  no  signal 
was  present.  Q*»*N  would  stand  for  a  particular  signal  model  in  the  presence 
of  noise  and  P(QCVk)  would  be  the  fraction  of  time  that  this  signal  is  present. 

A  statement  of  the  a  posteriori  probability  using  Bayes  theorem 
for  noise  alone  is 


PCQulXO 


P 

"Tair 


P(Qn') 


and  for  signal  plus  noise  the  a  posteriori  probability  is 

p  (  qC.n  \  xj  =  pu9- ' frvO  P(Ql^) 

1  P  Ui) 


From  Equation  II-  1 ,  one  has 
l 


PvXjQw)-  i 


and 


-p C-  i 


(II-  2) 


(II- 3) 


(H-4) 


(II-  5 ) 


with 
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•  2  r.  *•«  p  H  J  >■•} . 


(II- 6 ) 


From  this,  we  have  the  a  posteriori  probabilities, 


PCQuU^j  = 


P^k) 

I 

Lm) 


ex  p  [‘  i 

**  PUUJ 


C0«T'x4 


and 


ClT7)Wj-  PtXi) 


( II  -  7 ) 


(II-8) 


The  denominator  of  Equations  II- 7  and  II-8,  which  is  propor¬ 
tional  to  PCX;)  ,  normalizes  PCCMXi)  and  the  PC \Xi)  30  that  their 
sum  is  one.  For  Gaussian  distribution,  the  Bayes  estimates  of  11-7  and  II-8 
are  optimum  processes  since  the  overall  probability  of  error  is  minimized. 

Equations  II-7  and  II-8  are  the  basic  equations  for  their  multi- 
possibility  probabilistic  processor  using  Gaussian  distributions.  There  are, 
however,  two  other  essentially  equivalent  forms  of  Equations  TI-7  and  II-8. 
One  is  the  probability  ratio  form, 


pcj&ixi)  =.  itw «.pHx’rCGr^ -  o:i .  (II-9) 

Equation  II-9  gives  the  ratio  of  the  a  posteriori  probabilities 
of  P(Q«£U,)  and  P(Qn).  This  form  avoids  computing  the  normalization  factor, 
p(Xi)  used  in  Equations  II-7  and  II-8.  However,  this  normalization  factor 
can  beJ found  from  the  probability  ratios  since  the  sum  of  the  probabilities  is  one. 


The  second  equivalent  form  is  the  log  ratio, 
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9  pcq«Ui)  %jx^. 


(II- 10) 


This  form  shows  clearly  that  if 


then 

P(QwUx)>  PCQj^iU) 

(n-  ii) 


i.  e.  ,  the  data  is  more  likely  noise. 

Thus  the  problem  is  to  determine  which  is  more  likely,  Qn  or 
.  This  is  done  by  calculating  the  test  statistic, 


Xi  T  L  Q%4Vi  “  Owl  X* 


(11-12) 


from  the  datajX^  /or  each  plane  wave  signal  model  ■  The  number 

obtained  is  compared  with  the  threshold  value 


C  KV\  ~ 


z 


lo9c. 


lQ*/*PlQrtM\ 
WUl*  P CQn) 


(11-13) 


which  is  a  scalar  constant  for  a  particular  data  set  when  all  signal  models  are 
assumed  to  be  equally  likely.  Since  the  plane  wave  signal  models  have  cross¬ 
power  matrices, ,  which  are  functions  of  the  vector  wavenumber,  the 
test  statistic  of  Equation  II- 10  can  be  contoured  as  a  function  of  wavenumber.  In 
the  resulting  contoured  plot,  the  smaller  numbers  indicate  the  greater  probability 
that  a  plane  wave  signal  is  present. 


Sixth  Quarterly  Report 
Contract  AF  33(657)- 16678 


-14- 


1 1  January  1 968 


The  use  of  the  likelihood  ratio  test  statistic  of  Equation  II- 11 
is  favored  over  direct  use  of  the  a  posteriori  probabilities  of  Equation  II- 7 
and  II- 8  because  of  the  difficulty  in  establishing  the  a  priori  probabilities 
and  PtQ„)  .  Equation  II- 11  has  a  threshold  term  containing  the 
a  priori  probabilities  but  the  size  of  the  threshold  i3  not  of  paramount  im¬ 
portance.  The  relative  magnitude  of  Term  11-12  for  different  data  sets,  , 
or  different  signals  in  the  same  data  set,  contains  significant  information  as 
to  probable  event  location. 

B.  PROBABILISTIC  PROCESSING:  LARGE  ARRAY 

For  convenience,  the  likelihood  ratio  test  statistic  is  duplicated 

below, 


>  Z  'o9e 


lQJVt  PlQs~») 

\Q£v\'lL 


(HI  •  l) 


and  will  be  referred  to  as  the  probabilistic  processor. 

1.  Noise  Crosspower  Matrix  Estimation-Uncor related  Gaussian 
Noise  Field 

The  implementation  of  Equation  III-  1  becomes  simplified  when 
the  ambient  seismic  noise  is  assumed  to  be  uncorrelated  between  channels  and  to 
be  Gaussian.  For  a  large  array  such  as  LASA,  this  noise  field  assumption 
has  been  experimentally  shown  to  be  valid.  ^ 

For  the  above  mentioned  noise  field,  the  crosspower  matrices 
and  are  simplified  to 

-  Lvv'T-*  "^3  (m-2) 


-  t  *  xl 


and 


CQ«1 


(III- 3) 
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where 


f  JL  /l, 

X  t. 


With 


vector  wave  number 


f  n 

^  n  array  coordinate  corresponding 
to  the  data  channel 


W  = 


identity  matrix 

scalar  multiplier  equal  to  the  magnitude  of  the 
power  of  the  uncorrelated  noise 


The  superscript,  nq  ,  has  been  dropped  from  the  signal  plus  noise  crosspower 
matrix  since  the  plane  wave  signal  will  hereafter  be  implicitly  defined  as  a 
function  of  vector  wavenumber, 

The  inverse  crosspower  matrices  become 

CQww]  -  ijLi  “  (tv  +  V  V  * T  n 


(HI-4) 


ZQ-Y'  -  -b  Cl! 

Substitution  of  Equations  III-4  and  III- 5  into  the  left  side  of 
EquationfII-1  results  in  the  probabilistic  processor  assuming  the  form 


(III- 5) 


The  formula  used  is 

Ca^XY'T1  =  H-RX  (i+X'RY  )'V"R 

where  X  and  V  are  row  matrices  (vectors)  and  R  =  P\  Both  A  and 
(A+  XY**)  are  assumed  to  be  nonsingular. 


Sixth  Quarterly  Report 
Contract  AF  33(657)- 16678 


-16- 


11  January  1968 


vr  tQsi, -Q-  IXi  =  -  wTcWj  ( » 'XSHMS)  '  ? 


2  1©9 


1  P(Qs-*oj 

P  COo) 


(III- 6) 


After  some  algebraic  manipulation,  Equation  III  -  6  becomes 


t  Q*o  1 x  PCQv»vj) 
c\QvwlVlPCQM) 


I 


\n[w  +  \] 


^vj  (vTz')(i/Trs) 


7  o 


(III- 7) 


Equation  11-11,  of  the  previous  section,  states  that  if  this  inequality  is  true, 
then  more  probably  only  noise  is  present. 

The  calculation  on  the  left  side  of  inequality  is  conveniently  displayed 
as  a  function  of  discrete  vector  wavenumbers  in  a  format  similar  to  that  of  a 
frequency  wavenumber  spectra.  A  plane  wave  signal  would  appear  as  a  low 
magnitude  area  in  the  wavenumber  plot. 

2.  Equivalence  to  High  Resolution  Frequency- Wavenumber  Spectra  and 
Time  Domain  Beam  Steer 

High  resolution  frequency-wavenumber  spectra  are  computed  from 
a  frequency  domain  multichannel  filter.  In  the  filter  synthesis,  the  signal 
crosspower  matrix  is  assumed  to  be  uncorrelated  Gaussian  noise  and  the  noise 
crosspower  matrix  is  measured  data. 

Multichannel  filter  design  at  a  single  frequency  is 


=  r, 


L10I 


XiXSTl  fV 


(IH-8) 
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channel  is  reference 


where  K,  wO  ,  Xi  and  I  have  the  definitions  given  in  previous  sections.  The 
subscript,  n  ,  (lift  5«)  specifies  the  reference  channel  used  in  the  filter 
design. 


is 


The  solution  to  Equation  III- 8  using  the  n 


th 


channel  as  reference 


Fn  [T  ~{w  +  xS*Xl)  Xi  X^T1  rn 


(III- 9) 


A  high  resolution  frequency-wavenumber  spectra  may  be 
computed  from  Equation  III- 9  with  a  single  sensor  as  reference,  or 


V  r  Fn*  FS  1/  ' 


(III- 10) 


V  remains  unchanged  from  the  earlier  definition. 

Summation  of  the  responses  from  using  each  channel  as  re¬ 
ference  results  in 


\/rl2FSF„T J  ]/ 


(III- 11) 


Let  $  r\  ~  Pr)  m 
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Then  Equation  III- 12  becomes 


2  KrF„'F zv*  --  VTC  2  tS/fSlv' 

"w  t  yJrXL  ^  V 


(III- 12) 


Hcwever,  when  all  sensors  are  used  as  reference 

2  01*  =  I 

VI 

and  Equation  III- 12  simplifies  to 


(III- 13) 


(III- 14) 


Comparison  of  Terms  III-7  and  III- 14  reveals  that  the 
probabilistic  processing  equation  and  the  high  resolution  frequency-wavenumber 
spectra  (when  every  sensor  is  used  as  reference)  compute  a  test  statistic 
of  the  same  form.  That  is,  both  test  statistics  consist  of  a  constant  minu, 
amplitude  scaled  frequency  domain  beam  steer. 
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Since  the  presently  used  time  domain  beam  steers  are 
steered  for  a  fixed  velocity,  their  outputs  would  be  essentially  equivalent 
to  summing  over  like  velocities  the  probability  pictures  (high  resolution 
spectra  or  probalistic  processing  plot)  over  all  frequencies  in  the  bandwidth 
of  the  time  domain  processor. 

As  shown  above  for  a  random  noise  field,  probalistic  processing 
is  equivalent  to  beam  steering  or  the  high  resblution  spectra  using  all  sensors 
as  reference. 

C.  PROBABILISTIC  PROCESSING:  SMALL  ARRAY 

As  mentioned  in  an  earlier  section,  the  seismic  noise  field 
may  be  assumed  to  be  Gaussian  and  uncorrelated  between  sensors  for  larg 
arrays  such  as  LASA.  However,  as  the  array  aperture  is  decreased,  the  noise 
field  between  sensors  becomes  correlated.  The  simple  processor  described  in 
Section  III  is  no  longer  optimum. 

One  wo  uld  expect  signal  detection  capability  to  improve  for  small 
arrays  if  a  better  estimate  of  the  noise  crosspower  matrix  were  available. 

Since  the  noise  field  is  non-time-stationary,  the  estimate  should  vary  as  a 
function  of  the  prevailing  noise  conditions. 

1.  Noise  Crosspower  Matrix  Estimation  -  Measured  Data 

At  a  single  frequency,  a  time  varying  noise  field  can  be 
characterized  by  a  crosspower  matrix  computed  using  an  exponentially  weighted 
recursive  least  squares  scheme.  ®  Noise  croespower  matrix  updating  assumes 
the  form 


-  ^  [QutJ 


(IV- 1) 


where  t±t 

r 


time  interval  between  data  samples 
time  constant  of  data  decay 

interation  index  referring  to  a  particular  time 


Sixth  Quarterly  Report 
Contract  AF  33(657)- 16678 


-20- 


1 1  January  1968 


The  rate  at  which  past  estimates  of  the  noise  crosspower 
matrix  are  decayed  is  a  function  of  .  A  rapidly  varying  noise  field  would 
necessarily  require  a  more  rapid  taper  for  successful  tracking. 

Since  noise  field  measurements  in  a  signal  interval  are  biased 
by  the  presence  of  the  signal,  these  measurements  should  not  be  included 
in  the  noise  crosspower  matrix  estimation  of  Equation  IV-1.  A  decision 
criteria,  for  including  or  not  including  measurements,  based  on  amplitude  scan 
of  the  probabilistic  plot  could  be  used.  When  a  measurement  is  not  used  in  the 
crosspower  matrix  estimation,  the  matrix  LQ*+>jiw  will  remain  unchanged  from 

2.  Adaptive  Probabilistic  Processor 

Using  the  noise  field  estimate  of  Equation  IV-1,  the  signal  plus 
noise  crosspower  matrix  with  a  plane  wave  signal  model  becomes 


=-  +-\1Q  *0;,4i  ~\ 


(IV-2) 


- 1 

I 


Substitution  of  Equations  IV-2  and  IV-1  into  Equation  III- 1  produces 
an  adaptive  probabilistic  processor  which  is  optimum  when  the  noise  field 
As  non-time-s  <'tionary  and  correlated  between  channels.  The  processor  test 
statistic  is 


2  I03, 


!QjVa  PlQs-n) 

\QWhPCQ«) 


(IV- 3) 


1  +  i/'rC  Q 


i  '  '■ 


?  0 


As  is  evident,  the  inclusion  of  a  non-random  noise  field  complicates  the  com 
putation  of  the  test  statistic. 
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D.  PROBABILISTIC  PROCESSOR  PROGRAM  IMPLEMENTATION 

Wavenumber  displays  of  the  test  statistic  produced  by  a 
probalistic  processor  can  be  implemented  through  a  series  of  programs  re¬ 
presented  bv  the  flow  chart  shown  below. 


Individual  programs  of  this  series  can  be  combined  into  a 
single  program  which  will  take  multichannel  time  series  data  as  input  and 
sequentially  output  wavenumber  likelihood  ratio  plots. 
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E.  STUDY  OF  LOCATION  SCHEME  STATISTICS 

A  comparison  of  the  various  frequency-wavenumber  (f-k) 
location  schemes  in  several  signal  and  noise  environments  is  being  conducted. 
The  objective  of  this  study  is  to  determine  for  each  scheme  the  probability  of 
correctly  locating  plane  wave  events  in  both  correlated  and  uncorrelated  noise 
fields  at  varied  signal-to-noise  ratios.  With  this  information,  it  will  be  possible 
to  choose  the  best  scheme  for  locating  teleseismic  events  in  wavenumber  space. 

The  four  location  schemes  being  compared  are 

•  Convention  wavenumber  spectra 

vyX*V 


Probabilistic  processing  for  single  plane  wave 
signal 


rm  4v\rry»  x'k-'x- 


C  i  +  V/'KJ'Vj 


Two  forms  of  high  resolution  wavenumber  spectra 
HR  1 


V‘lN  +  xx"]'V  =  V  VV  - 


V'N'xX  *N~y 

(/  +  xv  u-'x) 


and  HR  2 


v'WtXxTV  -vW-  cy/* ~*V "fx&V 


hi  is  noise  covariance  matrix 
X  is  vector  of  data  transforms 
V  is  unit  plane  wave  transform  for  various 
signal  models 
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Initially  an  array  geometry  equivalent  to  the  subarray  locations 
(center  seismometer)  on  the  A,  B,  C,  and  D  rings  of  the  LASA  was  used  and 
the  signal  transforms  were  obtained  from  this  geometry.  The  noise  covariance 
matrix  N  becomes  a  constant  time  the  identity  matrix  pX  since  the  noise  is 
largely  incoherent  for  arrays  of  this  size.  6  in  this  situation  the  above  schemes 
reduce  to 


•  Conventional  wavenumber  spectra 

m'XX’V 

•  Probabilistic  processing 

\jv  y  _  X 

x  *  TTT~vrvJ 


\z*xx*v_ 

(  /  +  x'x) 


A  A 


•  HR  1 


i/V- 


HR  2 

v*v- 


L 


7-  4  1 

(  /  f  X*x)X  J 


V  'x  xV 


The  term  X¥X  is  the  sum  of  the  power  in  each  channel  at  the  frequency  being 
processed  and  ]/*!/  is  a  scalar  equal  to  the  number  of  channels.  Thus  the 
last  three  techniques  are  different  combinations  of  the  conventional  wavenumber 
spectra,  the  sum  of  the  power  in  each  channel,  and  the  number  of  channels. 

In  both  correlated  and  uncorrelated  noise  fields,  the  distributions 
of  each  of  the  four  techniques  cannot  be  expressed  in  an  analytically  closed 
form.  This  results  from  not  being  able  to  express  the  distribution  of  the  con¬ 
ventional  wavenumber  spectra  except  for  infinite  velocity  signal  models.  At 
infinite  velocity  the  distribution  of  a  conventional  wavenumber  spectra  is  the 
uniform  sum  of  chi-square  variables.  At  other  velocities  the  distribution  is  a 
non-uniformly  weighted  sum  of  chi-square  variables  whose  distribution  is  not 
known  in  closed  form.  Various  approximations  to  these  distributions  do  exist 
and  can  be  used  if  sufficient  experimental  data  is  available  to  define  the 
approximation  parameters.  This  study  will  provide  data  from  which  appropriate 
approximation  can  be  determined.  Finding  this  approximation  is  not  a  major 
objective  of  the  study,  but  it  may  be  very  useful  in  analyzing  the  relative  merits 
of  the  location  scheme  being  studied. 
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The  statistics  are  obtained  by  synthesizing  100  independent  trans¬ 
form  vectors  consisting  of  100  noise  transform  vectors  added  to  a  plane  wave 
signal  transform  vector  at  the  desired  signal-to-noise  ratio.  For  the  large  array 
the  noise  vector  was  made  up  of  zero  mean,  unit  variance  random  numbers 
generated  by  a  digital  computer.  This  simulated  an  uncorrelated  noise  field. 

The  data  from  this  experiment  is  being  evaluated  to  determine  the  locatability 
of  events  sensed  by  a  large  array. 

Further  distribution  studies  will  be  conducted  using  array 
geometries  of  the  extended  E3  subarray  and  a  standard  LASA  subarray  to 
determine  the  relative  performance  of  the  various  location  schemes  in  partially 
correlated  noise  fields.  The  results  of  this  study  should  indicate  the  relative 
merits  of  each  approach  when  operating  on  data  sensed  by  small  arrays. 

F.  CONCLUSIONS 


In  view  of  the  above  facts,  it  appears  unlikely  that  the  high 
resolution  wavenumber  spectra  techniqu  will  offer  any  improvement  over 
the  conventional  beam  forming  in  processing  large  seismic  arrays  such  as 
LASA  It  is  qui'-e  possible  that  computational  savings  are  available  by  beam¬ 
forming  in  the  frequency-domain  rather  than  the  time-domain.  These  com¬ 
putational  savings  would  be  a  result  of  the  numerical  methods  utilized  in  im 
plementing  the  present  techniques  rather  than  the  advent  o,  new  processing 

techniques. 

Degradation  of  beam  forming  performance  due  to  travel-time 
anomalies  and  signal  waveform  distortion  across  LASA  can  be 
especially  when  frequency-domain  processing  is  being  considered  Small 
arrays  such  as  the  extended  LASA  subarray  (E3)  and  a  standard  LASA  sub- 
array  exhibit  smaller  travel-time  anomalies  and  waveform  distortion  than 
does  the  large  array;  but  the  seismic  noise  field  is  correlated  between  individual 
elements  (seismometers)  of  the  array.  Observations  of  the  subarray  seismic 
noise  fields  (short  period)  indicate  that  they  are  non-isotopic  and  non- stationary. 
Th-s  ^tuation  might  require  an  adaptive  scheme  to  track  the  noise  field  m  order 
to  achieve  the  maximum  available  performance  from  small  array  processing. 

We  believe  that  ,t  would  be  fruitful  at  this  time  to  direct  our 
effort  toward  comparing  the  locating  ability  of  an  < iptimum 

abilistic  processor  with  the  location  ability  of  a  large  array  (LASA  using  the 
beam  forming  technique  and  with  high  resolution  f-k  spectra  of  small  arrays. 
TheTe  comparisons  will  be  conducted  using  moderate  and  weak  teleseismic  event, 
recorded  at  the  Montana  LASA  when  the  E3  extended  subarray  was  in  operation. 
This  will  permit  the  comparison  of  the  small  subarray,  the  extended  subarray, 
and  the  large  array  using  the  same  event  for  each. 
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The  application  of  probabilistic  processing  to  arrays  with  correlated 
noise  fields  will  require  some  modification  to  the  present  programs,  and  will 
not  be  as  numerically  efficient  as  a  processor  for  uncorrelated  noise  fields. 

The  increased  processor  complexity  might  be  offset  by  the  capability  to  accurately 
locate  events  using  array  of  small  aperture,  and  the  resultant  savings  in  array 
construction.  Determination  of  the  relative  capability  of  large  and  small  arrays 
when  processed  in  an  optimum  fashion  in  the  frequency  domain  is  then  a  reason¬ 
able  and  useful  objective  of  the  present  study. 

G,  SUBARRAY  PROCESSING 

To  prepare  the  raw  seismometer  data  for  the  detection  and 
location  processing,  the  events  being  used  are  processed  on  a  subarray  or 
individual  seismometer  basis.  For  the  large  array  processing,  the  sum  of 
the  seismometers  on  the  1,  3,  4,  5,  and  6  rings  of  each  subarray  is  formed; 
and  this  sum  is  bandpass  filtered  using  a  19-point  digital  filter.  The  band¬ 
pass  filter  has  a  sharp  low  cut  below  C.  65Hz  and  a  moderate  cutoff  above  1.  5 
Hz  and  is  designed  mainly  to  reduce  "leakage"  in  the  Fourier  transforms 
due  to  the  i licroseismic  noise  peak  at  0.3  Hz.  This  filter  will  also  be  applied 
to  the  individual  seismometer  outputs  before  the  extended  E3  and  standard 
LASA  subarrays  are  processed. 

The  particular  choice  of  seismometers  to  be  summed  in  the 
subarray  partial  sum  resulted  from  studying  the  array  response  of  several 
subarray  sums  as  well  as  the  response  of  an  MCF  designed  to  pass  greater 
than  10  km/sec  signals  and  reject  4  km/sec  noise.  The  MCF  did  not  appear 
to  significantly  outperform  the  chosen  partial  sum  to  warrant  its  use.  Details 
of  this  study  will  be  contained  in  Large  Array  Signal  and  Noise  Analysis 
Special  Scientific  Report  No.  18,  SUBARRAY  PREPROCESSING  FOR 
DETECTION  AND  LOCATION. 


Very  truly  yours, 

TEXAS  INSTRUMENTS  INCORPORATED 

Frank  H.  Binder 
Program  Manager 

FHB:  se 
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